[image1BkgRem,backgroundRow, imRowNoise, pRow, muRow, a] = remove_background(image1,5);
display_image(toggle,image1BkgRem,[0 100],'Image with background removed from row');
display_image(toggle,backgroundRow, [200 800], 'Background with stars removed (After Row fits)');
display_image(toggle,imRowNoise,[0 10],'Noise (row fit)');
[image1BkgRem,background1, imColNoise] = remove_background(image1BkgRem');
image1BkgRem = image1BkgRem';
imColNoise = imColNoise';
background1 = background1'; %Background of the image i.e. without the stars.
display_image(toggle, image1BkgRem, [0 100], 'Image with background removed from row and column');
display_image(toggle,background1, [0 50], 'Background with stars removed (After Row and Col fits)');
display_image(toggle,imColNoise,[0 10],'Noise (col fit)');
plot_aer_stars(realstar.locationAzEl(:,1), realstar.locationAzEl(:,2),...
realstar.brightness*50, 'r', 0, 0, 0, 1, 0);
plot_aer_stars(dascstarCal.locationAzEl(:,1), dascstarCal.locationAzEl(:,2),...
dascstarCal.brightness*50, 'c',...
calPar(1), calPar(2), calPar(3), calPar(4), calPar(5), calPar(6));
legend('Star chart','From calibrated parameters');
[dasc.azCal, dasc.elCal] = calculate_new_AzEl(dasc.az,dasc.el,calPar);
resize_figure(h,300,300);
p(1) = plot_DASC_aer(image1(indx), dasc.azCal(indx), dasc.elCal(indx), 1024, dsign);
p(2) = plot_aer_stars(realstar.locationAzEl(:,1), realstar.locationAzEl(:,2),...
realstar.brightness*80, 'c', 0, 0, 0, dsign);
p(3) = plot_aer_stars(dascStarAz, dascStarEl, dascstarCal.brightness*50, 'r', 0, 0, 0, dsign);
plot_grid_aer([0, 90], 22.5, 'm');
plot_star_names(realstar, 8, 'w', 0, 0, 0, dsign);
legend([p(1), p(2), p(3)], 'Calibrated image','Star chart','Stars from calibrated grid');
accuracy = (median(min(pdist2([dascStarAz,dascStarEl],...
realstar.locationAzEl(1:length(dascStarAz),:)))));
text(-90,-90,['Accuracy : ~',...
num2str(accuracy,2),'°']);
function [ASINew, background, ASINoise, pRow, muRow, a] = remove_background(ASI, nPoly)
% Remove background along the row
rowIndx = 1:1:size(ASI,2);
colIndx = 1:1:size(ASI,1);
ASINew = zeros(size(ASI));
background = zeros(size(ASI));
ASINoise = zeros(size(ASI));
pRow = zeros([length(rowIndx),nPoly+1]);
lchord = sqrt(512.^2-(abs(512-i)).^2);
fitRangeIndx =512-round(lchord)+1:1:512+round(lchord);
x = colIndx(fitRangeIndx);
y = ASItemp(i,fitRangeIndx);
if length(fitRangeIndx)>10
pRow(i,:)=polyfit(x,y,nPoly);
a(i).fitRangeIndx = fitRangeIndx;
background(i,fitRangeIndx) = polyval(pRow(i,:),colIndx(fitRangeIndx));
brightStarIndx = ASItemp(i,:) > (background(i,:)+1*stdRow(i,:));
if ~isempty(find(brightStarIndx)>0)
ASItemp(i,brightStarIndx) = nan;
ASItemp(i,fitRangeIndx) = interp_nans(ASItemp(i,fitRangeIndx)')';
ASINew(i,fitRangeIndx) = ASI(i,fitRangeIndx)-background(i,fitRangeIndx);
ASINoise(i,fitRangeIndx) = ASItemp(i,fitRangeIndx)-background(i,fitRangeIndx);
function ASI = remove_noise_spikes(ASI, sigma_n)
for i=3:1:size(ASI,1) - 2
for j = 3:1:size(ASI,2) - 2
B1 = ASI(i,j-1); B2 = ASI(i,j-2);
C1 = ASI(i-1,j); C2 = ASI(i-2,j);
D1 = ASI(i,j+1); D2 = ASI(i,j+2);
E1 = ASI(i+1,j); E2 = ASI(i+2,j);
singleSpike = (A > 2*s)*(B1 < 2*s)*(C1 < 2*s)*(D1 < 2*s)*(E1 < 2*s);
doubleSpike(1) = (A > 3*s)*((B1>2*s)*(B2<2*s))*...
not((C1>2*s))*not((D1>2*s))*...
doubleSpike(2) = (A > 3*s)*not((B1>2*s))*...
((C1>2*s)*(C2<2*s))*not((D1>2*s))*...
doubleSpike(3) = (A > 3*s)*not((B1>2*s))*...
not((C1>2*s))*((D1>2*s)*(D2<2*s))*...
doubleSpike(4) = (A > 3*s)*not((B1>2*s))*...
not((C1>2*s))*not((D1>2*s))*...
ASI(i,j) = not(singleSpike).*A;
ASI(i,j) = not(sum(doubleSpike)).*A;
ASI(i,j-1) = not(doubleSpike(1)).*B1;
ASI(i-1,j) = not(doubleSpike(2)).*C1;
ASI(i,j+1) = not(doubleSpike(3)).*D1;
ASI(i+1,j) = not(doubleSpike(4)).*E1;
function [starImage] = faint_star_extracter(ASI, sigma_n)
error('mSz-MSz has to be even');
starImage = zeros(size(ASI));
multiWaitbar('Faint Star Extraction...',0);
for i=1:1:size(ASI,1) - MSz
for j = 1:1:size(ASI,2) - MSz
M = ASI(i:i+MSz-1,j:j+MSz-1);
in = M(1+dSz:dSz+mSz,1+dSz:dSz+mSz);
out(1+dSz:dSz+mSz,1+dSz:dSz+mSz)=0;
sig_n = sigma_n(i+dSz+ceil(mSz/2),j+dSz+ceil(mSz/2));
mu_out = nanmean(out(:));
if mu_in > 2*sig_n && mu_out < 2*sig_n
starImage(i:i+MSz-1,j:j+MSz-1) = M;
multiWaitbar('Faint Star Extraction...','Increment',id);
function [newstar, I] = filter_stars(imstar, elMax, astrometryFileStr)
% Function sorts stars according to brightness, and defines a FoV of
% imstar: A structure of star extracted from image
% elMax : Stars below this elevation will not be considered
% astrometryFileStr : A text file written in the format necessary
% for astrometry web utility: http://nova.astrometry.net/upload
% newstar: A structure with stars sorted according to intensity
astrometryFileStr = []; % File containing text that can uploaded to astrometry
% Restricting the field of view.
imLength = imstar.size(1);
LengthperElevation = imLength/180;
p.min = round(imLength/2) - round(LengthperElevation*FOV/2);
p.max = round(imLength/2) + round(LengthperElevation*FOV/2);
selectedStarIndx =(imstar.location(:,1)>p.min &...
imstar.location(:,1)<p.max & ...
imstar.location(:,2)>p.min & ...
imstar.location(:,2)<p.max);
newstar.location = imstar.location(selectedStarIndx,:);
newstar.brightness = imstar.brightness(selectedStarIndx)';
% Sorting the stars based on its magnitude.
[newstar, I] = sort_star(newstar);
if ~isempty(astrometryFileStr)
dlmwrite(astrometryFileStr,round(newstar.location));
function [newstar, I] = sort_star(newstar)
% Function that sorts stars according to brightness
[newstar.brightness, I] = sort(newstar.brightness,'descend');
newstar.brightness = newstar.brightness';
newstar.location = newstar.location(I,:);
function realstar = get_actual_stars(stars,elCutOff,magCutOff,dx,dy,drot, dsign)
starfilter=stars.vmag<magCutOff & stars.el>elCutOff;
[x,y] = get_aer_stars(stars.az(starfilter), stars.el(starfilter), dx, dy, drot, dsign);
realstar.location = [x, y]; %pixel location
realstar.brightness = stars.relIntensity(starfilter);
realstar.locationAzEl = [stars.az(starfilter), stars.el(starfilter)];
realstar.name = stars.name(starfilter);
[realstar, I] = sort_star(realstar); %s
realstar.locationAzEl = realstar.locationAzEl(I,:);
realstar.name = realstar.name(I);
function [dascstar,x,fval] = calibrate_stars(realstar,dascstar,azOld,elOld)
dloc = round(dascstar.location);
dascstar.brightness = dascstar.brightness./max(dascstar.brightness);
lindx = sub2ind(size(azOld),dloc(:,2),dloc(:,1));
dascstar.locationAzEl = [azOld(lindx), elOld(lindx)];
minElFilter = find(dascstar.locationAzEl(:,2)>=22.5);
dascstar.locationAzEl = dascstar.locationAzEl(minElFilter,:);
dascstar.brightness = dascstar.brightness(minElFilter);
dascstar.location = dascstar.location(minElFilter,:);
ndasc = length(dascstar.brightness);
starIndx = true(1,ndasc);
% dx, dy translation of center point
% drot - rotation along azimuth
% dr - is field of view (radius)
% k - radial distortion coefficient
lb = [-10, -10, -180, -10, -100 ];
ub = [+10, +10, +180, +10, +100 ];
options = optimoptions('ga','MaxGenerations',2000,'MaxStallGenerations',500);
[y(k,:),fval(k),exitflag] = ga(@starDistance,nvars,[],[],[],[],...
x = [y(minIndx,1), y(minIndx,2), y(minIndx,3), dsign(minIndx),...
y(minIndx,4), y(minIndx,5)];
% Removing planets and other objects not in the star database
[x2, y2] = get_aer_stars(dascstar.locationAzEl(:,1),...
dascstar.locationAzEl(:,2),x(1),x(2),x(3),x(4),x(5),x(6));
D2 = pdist2([x2, y2], realstar.location(1:ndasc,:));
unmatchedObjects = dmin2 > 0.5; %Objects do not match if greater than 1 deg off
nMatched = sum(~unmatchedObjects);
accuracy = median(dmin2);
if nMatched>10 && accuracy<1
starIndx = ~unmatchedObjects;
[y2, fval2, ~] = ga(@starDistance,nvars,[],[],[],[],...
x = [y2(1,1), y2(1,2), y2(1,3), dsign(minIndx), y2(1,4), y2(1,5)];
disp(['Number of stars available should be more than 10. N = ',...
num2str(nMatched),'. Redoing the calculation']);
disp(['Accuracy is not good enough, it is ~ ', num2str(accuracy,2),...
'. Redoing the calculation']);
[dascstar,x,fval] = calibrate_stars(realstar,dascstar,azOld,elOld);
function dmin=starDistance(x)
[x1,y1] = get_aer_stars(dascstar.locationAzEl(starIndx,1),...
dascstar.locationAzEl(starIndx,2),x(1),x(2),x(3),dsign(k),x(4));
dascstar.newxy = [x1,y1];
D = pdist2(dascstar.newxy,realstar.location(1:sum(starIndx),:));